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Abstract. The last years have witnessed a rapid development of three-dimensional models of Type la supernova explosions. 
Consequently, the next step is to evaluate these models under variation of the initial parameters and to compare them with 
observations. To calculate synthetic lightcurves and spectra from numerical models, it is mandatory to follow the evolution 
up to homologous expansion. We report on methods to achieve this in our current implementation of multi-dimensional Type 
la supernova explosion models. The novel scheme is thoroughly tested in two dimensions and a simple example of a three- 
dimensional simulation is presented. We discuss to what degree the assumption of homologous expansion is justified in these 
models. 
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1. Introduction 

Type la supernova (SN la) explosions are commonly at- 
tribut ed to th ermonuclear explosions of white dwarf (WD) 
stars llHovle & Fow ler. I960) in a binary system. The pre- 
cise scenario of thes e events is, however, controver sial (for 
a recent review see iHillebrandt & N iemeve/ '20001) and it 
seems well possible that different mechanisms contribute to 
the SN la class. The final decision must be made by compar- 
ing theoretical models with detailed observations. Substantial 
progress has been achieved on both sides during the past years. 

The rapid development of three-dim ensional models of ther- 
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monuclear supernova explosions (e.g. Hillebrandt et al., 200(r, 

IReinecke et al., 2002a b c; Gamezo et al., 2003; Calder et aL, 
I2OO4I) leads naturally to the question whether they are capable 
of reproducing observed features of SNe la. Systematic param- 
eter studies on the basis of thr ee-dimensional SN la explosi on 
models have become possible dRopke & Hillebrandtll2004ai) . 

For the models by iReinecke et alJ ll2002alhlch synthetic 
light curves have been calculated with very encouraging re- 
sults dSorokina & BUnnikovl l2003h . However, for a thorough 
comparison of SN la models with observation via synthetic 
light curves and spectra, it is inevitable to follow the explo- 
sio n simulations up to hornq logous expansion. The models 
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( 2002alhld) reach only to f = 1 .5 s and 
1 20031) stayed below t - 2s, although it 
seems likely that the former are more advanced in the explo- 
sion stage, since there the initial flame shape is superposed by 
stronger spatial perturbations which enhance the development 
of Rayleigh-Taylor like instabilities and turbulence. 



A further motivation for evolving the explosion models to 
later times is to explore the effects of slower turbulent combus- 
tion regimes (the so-called distributed burning). These apply to 
low fuel densities reached after f ~ 1 .5 s by expansion of the 
WD. The present work focuses on reaching the stage of homol- 
ogous expansion and the issue of late burning will be addressed 
in a separate study. 

In the following we report on modifications of the code 
of IReinecke et al.llll999lf2002al) that enable us to follow the 
explosion model for much longer times than previous mul- 
tidimensional simulations. These novel methods are tested 
in two-dimensional models. For these as well as for an ex- 
ample in three dimensions the approach to homologous ex- 
pansion is studied. Our simulations are based on the single 
degenerate Chandrasekhar ma ss deflagration scenario (see 
Hille brandt & N iemever. '200(f). However, the methods devel- 
oped here can easily be applied to other models, featuring for 
instance a delayed detonation. 



2. Homologous expansion in SN la explosions 

A few seconds after ignition of the flame SN la explosions 
are expected to approach self-similar (homologous) expansion, 
which is characterized by a fluid velocity proportional to the 
radius. 



v{r) oc r. 



(1) 



This implies that the relative positions of the fluid elements 
do not change anymore. Consider the equation of motion in 



2 



F. K. Rbpke: Following multi-D SN la explosion models to homologous expansion 



Lagrangian formulation, 
dv 

p— = -Vf -pVd), (2) 
dt 

with P, p, and denoting the pressure, the density and the 
gravitational potential, respectively. From this equation it fol- 
lows that no relative change of the velocity of the fluid ele- 
ments co-moving with the expansion is equivalent to a van- 
ishing pressure gradient and evanescent gravitational force. 
Hydrodynamical interaction ceases if VP - and free expan- 
sion is reached for VO 0. 

Another criterion for the approach to homologous expan- 
sion can be derived from the energy balance of the explosion 
process. All energy liberated in the explosion is generated by 
thermonuclear reactions. While only a tiny fraction of this en- 
ergy is radiated away giving rise to the observable luminous 
event, most of it is used to overcome the gravitational binding 
of the WD star and to accelerate the ejecta. Due to expansion 
the ejecta are diluted and cooled. Thus, approaching homolo- 
gous expansion, the gravitational energy and the internal en- 
ergy should become small compared to the kinetic energy. 

3. Numerical model and modifications 

3.1. Explosion model 

Our numerical model is based on the scheme proposed by 
iReinecke et al.l dTggQ.) and inclu des the improvements de- 
scribed bv lReinecke et alJ(l2002al) . 

The hydrodynamics is modeled by the Euler equations 
with species conversion and appropriate source terms to take 
into account nuclear reactions. The numerical solution of 
these eq uations is based on the Prometheus implementation 
jFrvxell et al.. 1989) of the Piec ewise Parabolic Method (PPM) 
bv lColella & Woodward! ill 9 84 . 

The equation of state (EoS) models the degenerate matter 
of the WD star taking into account an electron gas that is de- 
generate and relativistic to a variable degree, the completely 
ionized nuclei following the Maxwell-Boltzmann distribution, 
a photon gas and electron-positron pair creation. 

The appropriate description of the nuclear reactions would 
be given by a full reaction network. Since this is too costly 
to run concurrently with the explosion model, a simplified 
description is c hosen. Following the approach suggested by 
iReineck e et al.l ( l2002al) five species are included, namely a- 
particles, '^C, '*0, ^"^Mg as a representative of intermediate 
mass elements, and ^^Ni as a representative of iron group nu- 
clei. The composition of the WD prior to the explosion is as- 
sumed to be a mixture of '^C and '^O. At the initially high den- 
sities burning proceeds to nuclear statistical equilibrium (NSE) 
composed of a-particles and nickel. Depending on temperature 
and density in the ashes, the fraction of a-particles and nickel 
changes. Once the fuel density drops below 5.25 x 10^ g cm"^ 
due to the expansion of the WD, burning is assumed to termi- 
nate at intermediate mass elements and below 1 x 10^ gem"-' 
burning becomes very slow and is not followed anymore. This 
may not be a satisfactory approach when the evolution is simu- 
lated to later times since it is possible that even the much slower 



flame contributes to the production of intermediate mass el- 
ements. However, this will not affect the production of iron 
group elements and may therefore be neglected in order to ob- 
tain a first-order estimate of the produced nickel mass and the 
total energy. But for synthetic spectra it may have a significant 
impact. We will address this issue in a forthcoming study and 
ignore it in the following, focusing on the hydrodynamical evo- 
lution of the models rather than on a detailed nucleosynthetic 
description. 

Since the range of relevant length scales in a SN la explo- 
sion (from the radius of the WD of ~ 1000 km to the width of 
the flame which is less than a centimeter) is by far too large to 
be resolved in numerial simulations, the flame evolution has to 
be modeled by an appropriate effective scheme. This is accom- 
plished by treating the flame as a discontinuity between fuel 
and ashes. We will call this model representation of the flame 
"flame front" in the following. 

In the so-called flamelet regime of turbulent combustion, 
which applies to large parts of thermonuclear SN explosions 
(Niemever& Wooslev, 1997), the flame is wrinkled by inter- 
action with turbulent velocity fluctuations. These stem from 
a turbulent cascade, which is generated by large scale buoy- 
ancy (Rayleigh-Taylor) instabilities. The interaction with tur- 
bulence increases the surface of the flame front and accelerates 
its propagation. Turbulence dominates flame propagation down 
to the Gibson scale, which is far below the numerical resolu- 
tion. Below the Gibson scale flame propagation proceeds stably 
in the cellular burning regime and therefore these scales have 
httle effect on the macrosco pic flame evolution (Ropke e t all 
l2003HRoDke etail l2004i^tf). 

According to these considerations the discontinuity de- 
scription not only averages the finite internal flame structrue 
but additionally defines a mean flame location neglecting the 
wrinkling of the flame at scales below the resolution of the 
computational grid. In this way it represents a "flame brush" 
rather than a smooth flame. To determine the turbulent burning 
ve locitv of the flame front the sub -grid scale model described 
byE lemever &Hillebrandtl ill 9951) is applied. 

The propagation of the flame f ront is implemented fol- 
lowing the the level set technique bv lOsher"& Sethian (.198^ . 
This technique was ex tended to deflagration flames by 
ISmihanovski et all ( 1 1 99"^ and a scheme (the so-called "pas- 
sive implementati on") applicable to SN la explosion models 
was developed bv lReinecke et alJlll999l) . The basic idea is to 
associate the flame front F with the zero level set of a signed 
distance function G 

r:={r|G(r) = 0), |VG| = 1. 

Its temporal evolution is prescribed by the so-called G-equation 



The velocity of the front motion is given by 

Di -V + Sill (3) 

for the passive implementation with v and n denoting the fluid 
velocity and normal vector to the flame front, respectively. The 
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turbulent burning velocity St = V^isgs is determined by the 
subgrid-scale energy q'sgs provided by the sub-grid scale turbu- 
lence model. With this approach the flame front is described 
as a sharp discontinuity and topo logica l changes are handled 
without problems (Reinecke et allfl999l) . 

3.2. Co-expanding computational grid 

To evolve numerical models to homologous expansion it is 
necessary to capture significantly growing length scales in the 
computational domain. One approach could be to start with a 
wide-spaced static computational grid and use adaptive mesh 
refinement to resolve the relevant parts of the domain. This ap- 
proach will, however, become increasingly expensive with evo- 
lution time, since the developing turbulence will fill large parts 
of the domain in which the grid has to be refined. In explicit 
schemes the maximum tim e evolution for each hydro step is 
given by the CFL criterion dCourant et all[l928i) . It states that 
to ensure stability of the hydrodynamics solver the time step 
needs to be chosen smaller than the time the fastest possible 
wave would need to cross a computational grid cell. Thus re- 
fining the grid drastically increases the number of hydro steps 
necessary to evolve the model. The level set method together 
with the sub-grid scale model provide a high accuracy even at 
relatively coarse computational grids. Therefore we do not ap- 
ply adaptive mesh refin ement. 

The simulations by iReinecke et alJ (l2002alblcl) were per- 
formed on a grid that was equally spaced in the center of the 
WD. To capture at least parts of the expansion of the WD, the 
width of the computational cells was exponentially increased 
in the outer regions resulting in a highly non-uniform grid. In 
the following we will refer to this implementation as the static 
grid approach. An advantage of this approach is that the reso- 
lution of the flame front is relatively fine as long as it stays in 
the uniform inner part of the grid. If this is true for the stage 
of the highest energy generation rate, numerical convergence 
of the models can be reached already for a grid with 256 cells 
per spatial dimension (Reinecke et al., 2002c). The static grid 
approach, however, has a number of disadvantages: 

- With the static grid it is very expensive to follow the explo- 
sion model to homologous expansion. 

- The non-uniformity of the grid is an obstacle to the 
implementation of the level set method. The neces- 
sary re-in itialization of the signed distance function G 
jSussman et al.. 1994) is difficult in the elongated grid cells 
in the outer regions. 

- Some more advanced sub-grid scale models tested in sim- 
ulations similar to that presented below (Schmidt et al., in 
preparation) require a uniform computational grid. 

- In the enlarged outer computational cells the flame front 
representation is very coarse. This problem is less severe 
in models that are confined to only one spatial octant 
(which was the case in all previous simulations), because 
the buoyancy-induced burning bubbles developed preferen- 
tially along the axes. This is not true if this artificial symme- 
try constraint is abolished and simulations of the full WD 
star are performed (kRopke & HillebrandL .20041}.) . 



Therefore we implemented a computational grid that co- 
expands with the explosion. Here the grid spacing is allowed to 
vary and the resulting discretization of the equations of hydro- 
dynamics is neither Eulerian nor Lagrangian. With this scheme 
it is possible to track the expansion of the WD during the SN 
la explosion with a uniform computational grid. This approach 
will be denoted as co-expanding grid in the following. It over- 
comes the drawbacks of the static grid approach. 

However, due to the uniformity and the expansion of the 
grid, the flame front will be less resolved at the stage of 
maximum energy generation if one starts out with the same 
grid spacing as in the inner parts of the static grid approach. 
Therefore it has to be tested which resolution is necessary in 
the new implementation to reach numerical convergence. 

Our approach is bas ed on the "moving grid" technique' 
by IWinkler et al.l ( Il984l) . This technique was already imple- 
mented in the original Prometheus code (Frvxell et all ll989l) . 
We will give a brief outline of the d eriv ation of the underlyin g 
equations following Mulleil ill 994 and IWinkler etal] lll984l) . 
Let (d/dt) denote the Eulerian derivative taken with respect 
to fixed coordinates in the laboratory frame and (D/Df) the 
Lagrangian derivative taken with respect to a definite fluid el- 
ement. The derivative with respect to fixed values of the mov- 
ing grid ("moving grid derivative") will be designated (d/dt). 
Analogous to the fluid velocity 

D 

Dt 

where r^. = r^iro, t) specifies the position of a definite fluid 
element, we define the grid velocity 



^grid — ^^''grid- 

Here, rgrid is the position of a definite set of grid coordinates. 
The relative velocity between fluid and moving grid is then 
given by 

I'rel - V- Dgi-id. 

Lagrangian and Eulerian coordinates are special cases of the 
moving grid for v„nd = v and Dgrfd - 0, respectively. The 
Eulerian and Lagrangian time derivatives of the density a of 
an extensive quantity are related by 

— a = — a + v\a, 
Df dt 

with V taken with respect to Eulerian coordinates. Similarly, 

—a = —a + Dgrid Vfl 
dt dt 

relates the moving grid derivative to the Eulerian derivative. 

Introducing the Jacobian Jf of the transformation between 
coordinates defining the initial volume of a fluid element dV^^j^ 
and the volume dVfiuid = •^f'^^nmd same fluid element at later 



' Originally, Wink ler et alJ il984b called this technique adaptive 
mesh. To avoid confusion with adaptive mesh refinemen t meth ods, we 
will denote it moving grid technique, following MuUea 
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times it is straightforward to derive the Euler expansion for- 
mula (e.g. Muller 1994 ) 



Dln/f 
Df 



Vi), 



which leads to the Reynolds transport theorem 

^ f (fldyfluid)= r 1^ + V-(i;fl)|dyfluid. 

Analogously one obtains the moving grid expansion formula 
dln7 



df 



where J now denotes the Jacobian of the transformation be- 
tween a moving grid volume dVgrid - -/fdVg^jj and its initial 
volume dy" . ,. 

grid 

This leads to the moving grid transport theorem (see e.g. 
lMulleJl9 94) 



+ V-(fl«) dVgrid. (4) 



Inserting the general form of a balance equation in Eulerian 
form 

d 

—a - -V ■ iav) + s(a), 
dt 

where s(a) denotes the source term, into Eq. one obtains 
I s(a) dVgrid 

I a dVgrid +1 V ■ \a(v - Vgrid)] dVj 



df 



— a dy„rid + 



grid 



df 



dS. 



With the appropriate quantities a and source terms s(a) this 
defines the set of equations of hydrodynamics that is solved 
numerically. 

4. Simulation setup 

The two-dimensional simulations we report on, were carried 
out on one quadrant assuming axial symmetry and mirror sym- 
metry to the other half space. In the three-dimensional simu- 
lations, only one spatial octant was calculated assuming mir- 
ror symmetry to the other octants. All simulations with co- 
expanding grid started with a computational domain of [2.02 x 
10* cm]^ discretized on the respective number of computational 
cells. Reflecting boundary conditions were applied on all sides 
of the domain. 

The explosion simulation was started with a WD near the 
Chandrasekhar mass limit in hydrostatic equilibrium. Its initial 
temperature was set to To = 5x10^ K. For the composition a 
50% mixture of '^c and "'O was chosen. The central density 
Pc was set to 2.9 X 10^ gem"-*. With these initial values a WD 
model was constructed by integrating the equations of hydro- 
static equilibrium using the EoS described above (see lReineckd 
l200ll) . An important point is that due to the hmited range of 



the EoS this integration was stopped as soon as a density of 
lO^-'gcm"^ was reached. This value was kept constant in the 
outer regions ("pseudo-vacuum"), which thus cannot maintain 
hydrostatic equilibrium. The resulting fluid motions in these 
parts of the computational domain leave the physical evolu- 
tion of the explosion unaffected, since here the masses are very 
small. As will be discussed below, they can, however, lead to 
numerical difficulties. A WD of the mass of 2.797 x 10^^ g 
(1.406 Mq) resulted from the setup procedure. 

The flame was ignited near the center of the WD in a 
sphere which was superimposed by toroidal perturbations to 
accelerate the growth of Rayleig h-Taylor like ins t abilitie s. This 
corresponds to the c3 model o f i Reinecke etaP ll2002cl) . This 
initial condition is rather artificial. However, the exact way 
in whic h the flame formation proceeds is no t yet precisely 
known dOarcia-Senz & Wooslevl ( 1 19951) and IWooslev et alJ 
(2004 provide approaches), and it seems likely that it cannot 
be resolved in the kind of simulations presented in the follow- 
ing. Therefore one has to find an averaged initial flame pre- 
scription. For the exemplary cases, our simple model will suf- 
fice, but the initial condition of the burning front remains an 
open question of SN la explosion models. 

An issue that deserves some consideration is the prescrip- 
tion of the grid expansion velocity. The goal is to follow the 
evolution of the WD and thus the radial velocity of the edge of 
the star has to be tracked. However, in our setup as described 
above, it is problematic to define the edge of the WD. Neither 
the pressure nor the temperature vanish here. A certain low 
density is also not a good indicator of the WD's edge since 
the overall density drops with expansion and especially the low 
density of the ashes can become comparable to that "outside" 
of the star. We therefore followed the 1 .4 Mq mass shell. This 
was achieved by defining a belt that includes 0.006 Mq around 
this mass shell and determining the mean fluid velocity inside 
this belt. According to this mean velocity the grid coordinates 
were expanded self-similarly so that the 1 .4 Mq mass shell was 
located at a fixed position at our computational grid. By track- 
ing a very high mass shell one can prevent the star from reach- 
ing the domain boundaries even in case of strong asphericities. 

5. Two-dimensional simulations 

5.1. Comparison with static grid calculations 

In Fig.n]the flame front morphologies in a co-expanding grid 
calculation are compared to a static grid simulation at differ- 
ent epochs of the explosion process. Both simulations start out 
from the same flame shape and were carried out on an [256]^ 
cell domain. Until f - 0.3 s the evolution of both flame fronts 
is almost identical. This is triviaUy the case because the grid 
spacing is the same in the part of the domain that is occupied 
by the flame front. In the co-expanding grid simulation the grid 
spacing is changed according to the velocities around the mass 
shell of 1.4 Mq. This shell is reached by the first shock wave 
due to the prompt incineration of the material behind the ini- 
tial flame after f ^ 0.3 s. From that moment on the grid keeps 
expanding self-similarly. Therefore its spacing at f = 0.6 s is 
larger than that of the inner part of the static grid. The better 
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Fig. 1. Comparison of the flame front evolution in a simulation with co-expanding computational grid (solid) with a static grid 
simulation (dashed). The right and the upper axis ticks indicate the location of every fifth grid cell for the co-expanding and the 
static grid, respectively. 
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Fig. 2. Total energy in the simulations with co-expanding and 
static grid. 
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Fig. 3. Energy production rate in the simulations with co- 
expanding and static grid. 

resolution of the flame front in the static grid case leads to a 
slight deviation of the flame front shape evolution (cf. Fig.[0. 
However, at f = 0.9 s the flame front in the static grid simula- 
tion has already entered the part of the domain with exponen- 
tially increasing grid spacing and consequently its outer parts 
are less resolved than the flame front in the co-expanding grid 
simulation. Therefore the flame front in the latter case develops 
more features in the outer parts (see the plots corresponding to 
f = 1.2 s and f = 1 .5 s in Fig. 

Summarizing the results plotted in Fig.ID we find that the 
flame front in the co-expanding grid simulation develops less 
structure in the inner parts and more features in the outer parts 
than the corresponding flame front in the static grid simulation. 
However, the overall flame front evolution is very similar 

The total energies produced by the two models are plotted 
in Fig. 121 Both models produced too little energy to account for 
a SN la explosion. This is expected for two-dimensional sim- 



ulations with the chosen flame setup. The co-expanding grid 
simulation produced slightly less energy than the static grid 
case. This effect can again be attributed to the different evolu- 
tion of the resolution of the flame front. Figure|3]shows clearly 
that the energy production rate reaches a maximum at around 
0.7 s (The initial high value of the energy generation originates 
from the instantaneous incineration of the material behind the 
initial flame.). At this time, however, the grid spacing at the 
flame front location in the co-expanding case is already wider 
than the one in the static grid approach. 

From our comparison we conclude that the co-expanding 
grid scheme provides a reliable and robust model of the SN 
la deflagration phase. Small differences compared to the previ- 
ously used static grid approach can be attributed to the different 
resolutions of the flame front at different epochs of the explo- 
sion process. 

5.2. Numerical convergence 

Numerical convergence is a basis to judge on the credibility 
of simulations. In the scenario under consideration, the flame 
front is unstable and affected by turbulence on scales that reach 
far below the resolution that is achievable in multi-dimensional 
simulations. From this point of view, numerical convergence 
in the details of the flame morphology cannot be reached. 
However, in our models convergence in global quantities (like 
the total energy generation or the production of burning prod- 
ucts) is expected to result from the interplay of large-scale 
flame front features with the sub-grid scale model. Ideally, a 
lack of resolution of large-scale features in the flame front rep- 
resentation should be compensated by an increased turbulent 
flame propagation velocity determined by the sub-grid scale 
approach. Of course, a certain threshold of resolution will need 
to be exceeded to reach this regime in our numerical implemen- 
tation. 

Our final goal is to simulate SNe la in three spatial dimen- 
sions. However, currently it is computationally too expensive 
to perform three-dimensional convergence tests. We therefore 
take the approach of Reinecke et al. (2002c) and perform the 
resolution study in two dimensions. From this we conclude on 
the behavior in three dimensions. For our convergence tests we 
conducted a series of two-dimensional simulations with resolu- 
tions of [64]2, [128]2, [256]2, [512]^, [1024]2, and [2048]^ grid 
cells. The flame fronts and densities after 1.5 s (when burning 
has ceased) are visualized in Fig. 0] The model with the low- 
est resolution (Fig. 0^) shows a peculiar behavior In the bet- 
ter resolved simulations the large-scale features resemble each 
other to a certain degree. But also here the details differ The 
buoyancy instability of the flame front gives rise to a nonlinear 
evolution of the large-scale flame features. The particular real- 
ization of this is sensitive to small perturbations and thus the 
different numerical noise in the models may lead to variations 
in the large-scale features. Thus, even in numerically converged 
models a small scatter of the global characteristics is possible. 

For the static grid simulation lReinecke et al 1 ll2002cl) found 
numerical convergence in the energy production for grid spac- 
ings of the inner part of the domain below A.*: » 10^ cm. 
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Fig. 5. Resolution of the flame front in co-expanding simula- 
tions with (a) [64]2 ceUs, (b) [128]^ cells, (c) [256]^ cefls, (d) 
[512]2 cells, (e) [1024]^ ceUs, (f) [2048]^ cells, and in static 
grid simulations with (g) [256]^ cells and an inner grid spac- 
ing of 10^ cm, (h) [256]^ cells and an inner grid spacing of 
7.9 X 10^ cm. 



However, for the late stages, when the flame front enters the 
nonuniform part of the domain, the energy production failed to 
converge in that test. 

This points to an important issue in our simulations. Neither 
in the co-expanding nor in the static grid approach the resolu- 
tion of the flame front is constant. For the latter the flame front 
resolution is constant only in the early stages of the explosion 
but decreases rapidly once the flame enters the non-uniform 
part of the computational grid. In the co-expanding grid ap- 
proach, the resolution of the flame front decreases steadily. This 
is demonstrated in Fig. |5] where the resolution of the flame 
front is plotted against time. In this plot, graph g corresponds 
to a static grid simulation that has reached numerical conver- 
gence according to Reinecke et al. (2002c). It is obvious that 
only co-expanding grid simulations with [1024]^ or more grid 
cells resolve the flame front better at all times. It needs to be 
emphasized at this point, that for the co-expanding grid model 
the situation is further complicated by the fact that the grid res- 
olution is adapted to the expansion of the star Thus, more vig- 
orous explosions will decrease the numerical resolution faster 

The total energies produced in our simulations with differ- 
ent domain discretizations are plotted in Fig.|6l The runs with 
[64]^ and [128^] cells show a peculiar behavior The simula- 
tions with a resolution of [512]^ cells and better seem to be nu- 
merically converged. It is interesting to note that the model with 
[2048]^ cells releases slightly less energy than the [1024]^ cells 
model. As discussed above, this can be attributed to the slightly 
different large-scale realization of the ffame front which is ob- 
vious from a comparison of Figs.|^ and|^. The model with a 
resolution of [256]^ cells is converged in the early stages but 
the energy production falls behind the higher resolved ru ns for 
t < 0.7 s. This is similar to what was found bv .Reinecke et alJ 
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Fig. 6. Total energy in simulations with co-expanding grid for 
different numbers of computational grid cells. 




Fig. 7. Evolution of the chemical composition in simulations 
with co-expanding grid for different numbers of computational 
grid cells. 

(l2002d) and can be attributed to the decreasing resolution of 
the flame front. 

The evolution of the chemical composition in our simula- 
tions is determined by the propagation velocity of the flame 
front and the fuel density at the flame front and thus closely re- 
lated to the expansion of the WD, which in turn depends on the 
flame front propagation. The temporal evolution of the compo- 
sition is plotted in Fig.0 It is obvious that in the models with 
[64]^ and [128]^ cells the chemical evolution deviates consid- 
erably from the better resolved runs. Here an insufficiently re- 
solved flame front burns less material but also leads to a slower 
expansion of the star (cf. Fig.|3. For the other models the pro- 
duction of iron group elements can be regarded as converged 
but the [256]^ cells simulation shows less intermediate mass 
elements, which are produced in later stages of the explosion. 




Fig. 9. Evolution of the turbulent energy in simulations with 
co-expanding grid for different numbers of computational grid 
cells. 



Fig. 10. Deviation of the angular averaged velocity from ho- 
mologous expansion at different times in the two-dimensional 
simulation with [1024]^ cells. 



This makes the simulation with [256]^ cells the most im- 
portant to judge on numerical convergence. Both the evolu- 
tion of the total energy and the chemical composition in that 
model suggest that the resolution of the flame front becomes 
too coarse to be numerically resolved at f 0.7 s. According to 
Fig.|5]this corresponds to a resolution of ~15 km. 

Our study implies that for two-dimensional simulations 
this is given for models with more than [512]^ cells start- 
ing from a resolution of better than 3.95 km. A model with 
[256]^ (with (Ax)ini = 7.9km) cells will give the correct pro- 
duction of iron group elements but slightly underproduce in- 
termediate mass nuclei. Since three-dimensional simulations 
explode stronger than their two-dimensional counterparts (cf. 
iReineck e et al.lEo0 2a^. a somewhat better initial resolution is 
desirable here. It needs to be emphasized that the resolution 
is not only crucial for numerical conv ergence, but also for the 
implementation of the physical model. IReineck e et al. '("2002b") 
find the most energetic explosions for multi-spot ignition sce- 
narios. To accommodate an adequate number of initial flame 
spots, one needs, however a rather fine computational grid. 

5.3. Reaching homologous expansion 

Figure|8ldepicts the temporal evolution of our SN la simulation 
with [1024]^ cells up to f = 10.0 s. Note that flie color-coded 
density range varies in the snapshots. The solid contour rep- 
resents the zero level set of the G-function. This is associated 
with the flame front at early times, but after f ~ 2 s the density 
of the WD has become so low due to expansion that burning 
ceases in our description. It is, however, interesting to note that 
the G - level set is still a good indicator of the interface be- 
tween unburnt and burnt material at later times (see snapshot 
at f = 5.0 s and compare to the snapshot at t - 10.0 s where 
the G - contour was omitted). The reason for this some- 
what surprising behavior is that the turbulent energy in our 
simulation decreases drastically after about 2 s This is shown 



in Fig.Owhere the total sub-grid scale turbulent energy is plot- 
ted against time. One has to keep in mind, however, that ^sgs 
increases with the size of the grid cells, which even moderates 
the decrease of the turbulent energy plotted in Fig. |9] Hence 
later than t ^ 2s q^gs becomes so small that the fluid velocity 
V dominates in Eq. (|3} and consequently the G = isosurface 
(isoline) is advected passively with the resolved flow. 

Figure|9lagain demonstrates numerical convergence for our 
models with more than [256]^ cells. Here the criterion is not 
an exact match of the value of the turbulent sub-grid scale en- 
ergy but its temporal evolution. This is similar for the highly re- 
solved simulations while the less resolved runs diverge at early 
times. 

The morphology of the density structures does not signifi- 
cantly change after f ~ 5 s. This indicates that our simulation 
approaches homologous expansion. How consistent our simu- 
lation at late times really is with the assumption of homologous 
expansion has to be decided on the basis of the velocity profile 
and the amplitude of the gradients of pressure and the gravita- 
tional potential. 

This will be analyzed for our simulation with [1024]^ cells. 
Fig.^Jshows the deviation of the angular averaged velocities 
from relation Q at different times. The proportionality factor 
of that relation was obtained from a fit to the velocity profile 
at f = 10 s. The velocity deviations are plotted against the 
distance from the WD's center. This distance is expressed in 
widths of computational cells Ax to allow comparison; because 
of the uniform grid the physical distance scales directly with 
Ax. We note that the deviations from homologous expansion 
decrease drastically. Except from the region in the very center 
the deviation at f = 10 s stays below 5%. The static grid simu- 
lations ended att- 1.5 s. Here deviations still reach 10%. The 
maximal values and the mean values of the pressure gradients 
are given in Table[2 These quantities also decrease drastically 
with time. From t - 0.6 s (where the energy generation rate 
peaks) to f = 1.5 s they differ less than a magnitude while from 
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Fig. 11. Temporal evolution of the ratio of the kinetic energy 
to the sum of internal and gravitational energy for the two- 
dimensional model with [1024]^ cells (solid curve) and the 
three-dimensional model with [256]^ cells (dashed curve). 

Table 1. Pressure gradients and gradients of the gravitational 
potential in the two-dimensional [1024]^ cells simulation at dif- 
ferent times. 



t[s] 


max \Wp\ 


<|Vp|> 


max |V<1)| 




[dyn cm"'] 


[dyn cm"'] 


[cm s"^] 


0.6 


9.44 X 10' 


1.17 X 10' 


8.17 X 10' 


1.5 


6.85 X 10' 


2.95 X 10* 


4.05 X 10* 


5.0 


3.82 X 10"* 


8.74 X 10'' 


1.74 X 10' 


10.0 


1.27 X 10"* 


3.40 X 10" 


4.09 X 10** 



t - 1.5 s to f = 10.0 s they drop two magnitudes. The evolu- 
tion of the gradients of the gravitational potential (see Tabled 
shows the same trend. 

We conclude that our simulations at t a 10 s provide a fair 
approximation to homologous expansion. This is supported by 
Fig-[n(solid curve) showing that at f « 10 s the kinetic energy 
of our model is almost one order of magnitude larger than the 
sum of the internal and gravitational energies. 

Due to the co-expanding grid the computational costs to 
reach homologous expansion are moderate. The expanding grid 
spacing leads to an increase in the time steps satisfying the CFL 
criterion. For the Godunov-scheme used to solve the equations 
of hydrodynamics in our implementation this criterion states 
for the time step At 

At < Ax (max; + Cs,,})"' , 

where the index / runs over all computational cells and v and Cs 
denote the fluid velocity and the sound speed, respectively. The 
CFL criterion is necessary but not sufficient to ensure stability 
of the scheme. To guarantee stability we reduce the time step 
by a factor of 0.8 in our implementation. 

The evolution of the time steps taken in our simulation 
with [256]^ cells is plotted in Fig.^] Two points should be 
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Fig. 12. Time steps in the simulation with [256]^ cells. 

noted here. First, the time step is very low between ~0.30 s and 
~0.38 s. The reason for this is that the shock wave generated by 
the initial flame setup reaches the WD's surface and propagates 
into the surrounding pseudo-vacuum. Here high velocities are 
generated so that computational cells outside the star dominate 
in the CFL criterion. However, once the first shock wave has 
crossed the WD, it starts to expand and the computational grid 
co-expands accordingly. With the reflecting outer boundaries 
this leads to a damping of the motions in the cells filled by the 
pseudo-vacuum and the time step recovers to larger values. The 
second point to note here is that the increase of the time step 
due to the expansion makes it possible to reach from 2 s to 10 s 
with only ~10% more hydro steps than calculating up to 2 s. 

6. A three-dimensional simulation 

With the scheme described in the previous sections we per- 
formed a three-dimensional test simulation. The computa- 
tional domain was discretized into [256]^ cells starting with 
a grid spacing of 7.9 x 10^ cm and the flame setup coiTe- 
sponded to the three-dimensional version of the c5-model by 
iReinecke et al.l (k2002b,.c) . Snapshots from the simulations are 
shown in Fig. [21 The rendered isosurface coiTesponds to the 
zero level set of G. As discussed above it represents the flame 
front at early stages and still provides an approximate separa- 
tion between fuel and ashes after f ~ 2 s. 

The early evolution of the flame front is very similar to that 
found by Reinecke et al. ( 2002b c) . The initial axial symmetry 
is broken in the same way. The late stages show more structure 
due to the better resolution of the co-expanding grid in the outer 
parts of the domain. Our study in the previous section indicated 
that the simulation is close to homologous expansion already 
after 2 s. In the three-dimensional simulation, the morphology 
of the G - isosurface still changes. The outer bubbles grow 
hiding the inner structures. This may have impact on the light 
curves and the spectra derived from such models. From f = 5 s 
to f = 10 s the visible changes are only marginal. 




Fig. 13. Three-dimensional SN la simulation. The isosurface corresponds to G = and the color-coding indicates the value of 
the sub-grid scale energy providing a measure of the propagation velocity of the front. The scale gives an impression of the 
expansion and corresponds to the plane through the origin of the coordinate system. 
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Fig. 14. Comparison of the total energies of our three- 
dimensional simulation with co-expanding grid and the 
c3_3d_256 model (static grid) bv .Reinecke et al..c2002b.c.) . 
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Fig. 15. Deviation of the angular averaged velocity from ho- 
mologous expansion at different times in the three-dimensional 
simulation. 

The en ergetics of the m odel is compared to the c3_3d_256 
model^ bv lReinecke et alJ (|^002b c) in Fig.[2] As in the two- 
dimensional simulations the differences between the static grid 
and the co-expanding grid simulations can be attributed to the 
varying resolution of the flame front in the latter. 

Fig-Elshows the deviation from a self-similar velocity pro- 
file in our three-dimensional simulation. We note that it ap- 
proaches a relation consistent with Eq. ([0 faster than the two- 
dimensional simulations. This can at least partly be attributed 
to the more energetic explosion resulting in three dimensions, 
giving rise to higher expansion velocities. Figure [TTI corrobo- 

^ The data was kindly provided by M. Reinecke. In this model the 
grid spacing at the inner part of the computational domain was Ax = 
10^ cm. 



3.0 X lO--* 
2.5 X lO--* 

^ 2.0 X 10-^ 

's 







1 1 1 1 1 1 1 1 


lO.Os 








5.0 s 
1.5 s 






Mh- '^C/'-OA'. \ 
: II 1 ' V*, ■ * 












/"' 


■' ^ / V-., 1 \^ \ 






:' / 1 , 
/ / ■' 

/ ; 


1 \, \ \m 
1 \, 1 Y.\ 




j 

r 


/■7 


, , 1 , , 1^"" ^ 





3.0x10" 6.0x10* 9.0x10* 1.2x10' 

u,ad [cms"'] 

Fig. 16. Distribution of the species in radial velocity space. 
"Ni" and "Mg" denote the iron group elements and the inter- 
mediate mass elements, respectively. 

Table 2. Pressure gradients and gradients of the gravitational 
potential in the three-dimensional simulation at different times. 



r[s] 


max IVpl 


<|Vp|> 


max |V<D| 




[dyncm"'] 


[dyn cm"'] 


[cm s"^] 


0.6 


5.49 X 10' 


4.02 X 10* 


6.93 X 10' 


1.5 


5.57 X 10* 


1.50 X 10* 


3.16 X 10** 


5.0 


1.45 X 10^ 


4.74 X 10^ 


1.49 X 10' 


10.0 


5.70 X 10-* 


1.90 X 10^ 


3.59 X 10* 



rates this interpretation showing an accelerated conversion of 
the energy forms to kinetic energy in the three-dimensional 
model. The gradients of the pressure and the gravitational po- 
tential are given in Table |2] The maximum value of IVOI de- 
velops similar to the two-dimensional case while the mean 
of the pressure gradient amounts to roughly half the value 
obtained in the two-dimensional simulation. As in the two- 
dimensional model, the pressure gradient drops significantly 
between t - 1.5 s and f = 5.0 s. 

The distributions of chemical species in radial velocity 
space for different times are shown in Fig. ^] The changes 
between f = 5.0 s and f - 10.0 s are small. The profile is some- 
what shifted toward smaller velocities which corresponds to a 
deceleration of the ejecta. This is not surprising since the grav- 
itational potential is small, but not vanishing. In contrast, be- 
tween t - 1.5 s and t - 10.0 s we observe substantial changes 
in the profiles. All three species tend to pile up at lower veloci- 
ties diluting the high- velocity components. 

We conclude that our scheme is applicable to three- 
dimensional simulations. The expansion of the computational 
grid according to the tracked I.4M0 mass shell did not cause 
any difficulties. Homologous expansion in three-dimensional 
simulations will be reached somewhat quicker than in two- 
dimensional models. 
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7. Conclusions 

We analyzed the approach of multi-dimensional SN la mod- 
els to homologous expansion. To this end a novel scheme with 
co-expanding computational grid was implemented and tested. 
Taking into account the variable resolution resulting from this 
scheme, a comparison with previous simulations on a static 
computational grid showed good agreement in the flame evolu- 
tions. 

A convergence study validated our co-expanding grid ap- 
proach. Here, the variability of the numerical resolution leads 
to the conclusion that the grid spacing has to stay below ~15 km 
for numerically converged stages of the explosion. For two- 
dimensional simulations one has to start with a domain dis- 
cretization finer than [512]^ cells to resolve the whole explo- 
sion process, but a [256]^ cell models will already provide the 
correct production of iron group elements. In case of three- 
dimensional explosion models a somewhat finer initial resolu- 
tion is required because the explosions here are more vigorous 
resulting in a faster expansion. 

Applying our co-expanding grid scheme, it was possible 
to follow the explosion significantly longer than in any pre- 
vious multi-dimensional SN la simulation with only marginal 
additional computational expense. From the temporal evolu- 
tion of the velocity profiles, the energy conversion and the 
gradients of pressure and gravitational potential we conclude 
that homologous expansion is reached to a high accuracy at 
t ~ 10 s. Although the velocity profile is rather quickly in 
agreement with Eq. Q, the pressure gradient drops signif- 
icantly at later times. Therefore the flame morphology still 
changes until f ~ 5 s. Outer bubbles grow and hide inner parts 
of the ejecta. Between t - 1.5 s and t - 10 s we observe a 
significant redistribution of chemical species in velocity space. 
Thus, for the derivation of synthetic light curves and spectra 
from SN la models it is advisable to follow the evolution of the 
explosion up to f ~ 10 s. 

The novel numerical scheme presented here offers a variety 
of improvements of multi-dimensional SN la explosion mod- 
els. A uniform cartesian grid alleviates some numerical prob- 
lems (see Sect. ^ and moreover enables the implementation 
of more sophisticated sub-grid scale models (Schmidt et al., in 
preparation). Full star SN la explosion models, where the artif- 
ical symmetry to the spatial octants is a bolished, also gain ac- 
curacy from such a computational grid jRopke & Hillebrandtl 
l2004bV It is, however, also possible to combine the co- 
expanding grid approach with the highly non-uniform grid used 
in previous static-grid simulations. In this way a high resolution 
of the flame front (especially in the early stages) will be pos- 
sible with low computational expenses and initial flame con- 
ditions can be tested easily. We will address this question in a 
subsequent study. 

An open issue in the current implementation is the mod- 
ehng of flame propagation in late stages of the explosion. 
Here the finite width of the flame is not negligible anymore 
since turbulent eddies start to penetrate the preheat zone of the 
flame. The so-called flamelet assumption of turbulent combus- 
tion breaks down and the flame enters the thin reaction zones 
regime (for a discussion of the turbuelnt burning regimes see 



e.g. IPeteri fl995l and for ap plication to SN la explosion see 
Niemever & Kerstei'J Il997l) . It is, however, unclear whether 
burning in this regime contributes significantly to the explosion 
results. This issue will be addressed in a forthcoming study. 
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